AD-A094  771 
UNCLASSIFIED 


AIR  FORCE  INST  OF  TECH  WRI6HT-PATTERS0N  AFB  OH  SCHOO— ETC  F/®  ®/2 
COMPUTER  SIMULATION  OF  SOLAR  AXR  HEAT  I  NO  SYSTEMS  USING  ROCK  BED~ETC<U> 
DEC  80  D  B  FANT 

AFIT/6AE/AA/80D-4  NL 


AIR  UNIVERSITY  (ATC) 

AIR  FORCE  INSTITUTE  OF  TECHNOLOGY 


Wright-Pottorson  Air  Force  Base,  Ohio 


APPROVED  FOR 


23  JAw  j-j  o 

p::OLl':  LTtfZ  APR  190-17. 


!  'TAasA  <>  /  s  4V\°  ‘i  .  ■_ 


/  //  L 

LAUREL  A.  LAMPELA,  2Lt ,  USAF 
Deputy  Director,  Public  Affairs 


A:r  f c:cc-  L  '■ 

Wdgiil  i'  „UCi-„ 


c;  (ATC) 

i  hi  o,  oil  4o433 


[  Accession  for 

1 

I  N1IS  GRA&I 

; 

DTIC  TAB 

□ 

Unannounced 

n 

Justification _ 

i 

By.  . 

Distribution/  | 

Availability  Codes  ~j 
i  ail  and/or 

•List  StHcial  ! 


i 


COMPUTER  SIMULATION  OF  SOLAR 
>) I R  HEATING  SYSTEMS  USING 
ROCK  BED  THERMAL'  STORAGE  JJN ITS 


THESIS 


u  t  a.: 

ELECT  •  Vi, 

FEB  1  0  1981 


AFIT/GAE/ AA/80D- 4  - 


81  2 


09  094 


AFIT/GAI-/AA/S0D-4  ' 


computfr  simulation  of 

SOLAR  AIR  I  HATING  SYSTEMS 
USING  ROGK  Bl.D  TIIHRMAL  STORAGF.  UNITS 


THUS IS 


Presented  to  the  Faculty  of  tne  School  of  Engineering 
of  the  Air  Force  Institute  of  Tecimologv 
Air  Training  Command 
in  Partial  Fulfillment  of  the 
Requirements  for  the  Degree  of 
Master  of  Science 


by 

Daniel  B.  Fant 
2Lt  USAF 

Graduate  Ae ronaut  i ca  1  F.ng i nee r i ng 
December  1980 


Approved  for  public  release,  distribution  unlimited. 


Preface 


Rapid  interest  in  solar  energy  has  come  about  as  a  result  of  increasing 
costs  of  energy  from  conventional  fuels.  However,  due  to  the  intermittent 
and  diffuse  nature  of  solar  energy  some  unique  problems  have  also  arisen.  One 
in  particular  is  thermal  storage,  and  this  motivated  the  study  of  solar  air 
heating  systems  utilizing  rock  beds  as  thermal  storage  units.  A  computer  simu¬ 
lation  model  capable  of  estimating  the  long-term  performance  of  an  air  heating 
system  is  described  in  this  work,  along  with  complete  instructions  on  how  to 
operate  the  program.  Various  systems  can  be  simulated  by  simply  making  the 
appropriate  change  in  input  parameters. 

I  wish  to  thank  the  Air  Force  for  accepting  me  into  this  program  and  pro¬ 
viding  the  exceptional  staff  and  faculty  which  made  this  educational  experience 
most  worthwhile.  I  also  wish  to  thank  Mr.  Stan  Boyd  for  obtaining  research 
material  that  was  not  readily  available  from  the  AFIT  library. 

I  would  especially  like  to  thank  Dr.  James  E.  Hitchcock,  my  thesis  advisor, 
for  his  many  hours  of  assistance,  encouragement  and  suggestions  which  made  this 
thesis  possible.  But  most  of  all  I  would  like  to  thank  my  parents  for  their 
reassurance  and  support  which  in  so  many  ways  kept  me  going  and  made  me  strong. 


Daniel  B.  Fant 
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ABSTRACT 

This  thesis  is  concerned  with  the  analysis  and  design  of  solar  air 
heating  systems  utilising  rock  beds  as  thermal  storage  units.  A  computer  simu'ailo 
model  capable  of  estimating  the  response  of  both  the  solar  collector  an'5  tbe 
rock  bed  is  described.  r'ifrerential  equations  describing  the  rnofr  *>ed  were 
approximated  in  a  finite-difference  form  and  solved  numerien'Py  or  o 
computer,  ^he  temperature  of  hot’"1  tbe  solid  (reel'd  end  tHp  oRga  (air'  is 
determined  as  a  function  of  time  and  distance  along  ^he  bed.  ^e  Simula4"’ or 
required  both  charging  and  discharging  of  the  roc!'  bed  por  time-varying  ini ot 
fluid  tempertures.  The  numerical  method  used  to  solve  the  rock  hod  equations  prove 
to  be  stable  and  convergent  and  showed  satisfactory  agreement  in  comparison 
to  an  analytical  solution  for*  constant -inlet  air  temperatures. A  cost  analysis 
was  also  incorporated  within  this  program,  by  varying  the  collector  area  one 
could  determine  the  optimum  collector  sine  for  maximum  savings.  Pressure  drop 
relationships  for  flat-plate  collectors,  duct  work  and  packed  beds  were  used  to 
determine  operating  costs.  The  particular  air  system  tested  proved  to  be  cosi- 
effective  when  compared  with  natural  gas  fuel  costs  for  an  economic  term  of  on 
years . 
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COMPUTER  SIMULATION  OF  SOLAR 
AIR  HEATING  SYSTEMS  USING 
ROCK  BED  THERMAL  STORAGE  UNITS 


I  Introduction 


Background 

19 

The  incident  solar  energy  on  earth  is  approximately  1.2  x  10  Btu  per 
year,  while  the  world's  total  energy  consumpt ion  is  about  7  x  10^  Btu  per 
year  (Ref  17).  However,  solar  energy  is  a  variable  source  of  energy.  Here 
on  earth,  the  sun  shines  only  during  the  day;  and  on  cloudy  days  it  is  inter¬ 
mittent.  Therefore,  proper  storage  of  this  abundant,  cheap, and  pollution 
free  energy  is  necessary  for  its  efficient  utilization. 

The  idea  of  thermal  storage  goes  way  back  to  the  caveman  or  Roman  era 
when  they  brought  sun  heated  stones  or  bricks  to  bed  to  supply  them  with 
warmth  during  the  night.  The  Romans  also  dug  deep  cellars  in  the  side  of 
hills  which  they  filled  with  snow  during  the  winter  and  it  would  remain 
throughout  the  summer  to  be  used  when  needed.  These  above  happenings  can 
best  be  explained  by  the  fact  that  rocks  act  as  a  good  storage  material  (Ref  17). 

The  direct  mode  of  operation  of  air  heating  systems  is  to  transfer  heated 
air  directly  from  the  collector  to  the  heated  space.  When  more  energy  is  col¬ 
lected  than  needed,  it  is  stored  sensibly  in  a  storage  bed  by  blowing  heated 
air  from  the  collector  through  the  bed.  This  is  referred  to  as  the  bed  charging 
mode.  When  energy  is  not  available  from  the  collector  to  satisfy  the  space 
heating  loads,  sensible  energy  is  removed  from  the  rock  bed.  This  is  blown  as 
the  bed  discharging  mode.  Discharging  is  normally  accomplished  by  a  reversal 
of  the  air  flow,  thus,  causing  the  outlet  air  temperature  to  approach  the  tem¬ 
perature  at  the  hot  end  of  the  bed. 
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Purpose 

The  purpose  of  this  thesis  was  to  develop  a  computer  program  which  com¬ 
pletely  simulate*  the  three-modal  operation  of  a  solar  air  heating  system 
utilizing  rock  beds  for  thermal  storage.  This  involve*  a  numerical  method 
for  computing  the  temperatures  of  both  the  air  and  rock  as  functions  of  time 
and  distance  for  the  charging  and  discharging  modes  of  the  bed. 

With  a  proper  theoretical  analysis  of  rock  beds,  duct  work  and  collectors, 
a  simulation  model  would  provide  the  useful  purpose  of  estimating  the  performance 
of  a  solar  air  heater  much  more  quickly  and  cheaply  than  would  be  possible  by 
experimentation.  In  order  to  simulate  various  systems,  a  change  in  only  certain 
parameters  is  necessary. 

Approach 

Before  simulating  an  air  system  along  with  the  thermal  performance  of  a 
rock  bed  storage  unit  with  varying  inlet  air  temperatures,  one  must  first  develop 
a  scheme  for  solar  energy  collection.  Fortunately,  this  was  accomplished  by 
utilizing  the  previous  work  of  Captain  Prins  (Ref  16)  on  solar  water  heating 
systems.  His  long  term  analysis  of  flat-plate  collectors  was  used  for  deter¬ 
mining  the  useful  amount  of  solar  energy  collected  per  unit  time.  Here,  long 
term  is  defined  to  be  that  period  of  time  necessary  to  describe  the  performance 
of  a  heating  system  in  any  particular  climate,  normally  a  one  to  two  year  period 
is  sufficient  for  Jong  term  simulations  (Ref  12). 

The  amount  of  useful  energy  collected  is  the  difference  between  the  solar 
energy  absorbed  and  the  energy  lost  by  the  collector.  It  is  a  function  of  the 
collector  area  (AC),  the  collector  heat  removal  efficiency  (FR),  the  monthly 
average  instantaneous  total  radiation  (T.j,  ) ,  the  util  izabil  ity  (Iff)  and  the 
effective  transmittance-absorptancc  product  (m). 

Q  =  AC  x  FR  x  Tt  x  irr  X  ha)  (1) 
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FR  is  a  function  of  the  fluid  flowrate  and  the  absorber  plate  design. 
depends  on  the  location  and  orientation  of  the  collector,  also  on  the  hour 
of  day  and  the  time  of  year.  The  utilizability  is  a  function  of  variableness 
of  the  incident  solar  energy  and  the  loss  coefficient  of  the  collector,  which 
in  turn  is  a  function  of  the  number  of  cover  plates,  cover  material  and  amount 
of  insulation.  Lastly,  (toT)  also  depends  on  absorber  plate  properties,  the 
number  of  cover  plates  and  the  angle  at  which  solar  radiation  strikes  the 
collector. 

The  useful  energy  collected  in  equation  (1)  also  determines  the  temperature 
rise  of  the  collector  air. 


Cm  c  )  (T  -  T.  ) 
a  pa'  v  oc  ic  ' 


(2) 


where  rti  is  the  mass  flowrate  of  air  into  the  collector,  c  is  the  specific 
a  ’  pa  1 

heat  of  air  at  the  constant  pressure,  is  the  air  temperature  into  the 
collector  and  T  is  the  air  temperature  out  of  the  collector.  Therefore,  if 
Q  is  determined  from  equation  (1),  one  can  then  solve  for  Tqc  from  equation  (2)  . 


T  =  [Q/mc  )  !  +  T.  (3) 

oc  p'a  ic 

liquation  (3)  is  utilized  in  the  program  to  determine  the  varying  air  temperatures 
into  the  rock  bed  during  the  charging  mode. 

Flat-plate  collector  results  must  be  combined  with  a  rock  bed  performance 
analysis  to  completely  simulate  air  heating  systems.  This  involved  describing  a 
control  strategy  for  the  various  operational  modes  of  the  system;  along  with 
developing  an  adequate  space  heating  model  and  determining  auxiliary  energy 
requirements . 

A  slight  mod i fication  of  Prins'  cost  analysis  was  also  necessary  to  deal 


with  air  system  simulations. 


The  heart  of  this  thesis  remains  to  follow,  it  deals  with  the  theoretical 
heat  transfer  analysis  of  rock  beds  in  Chapter  II,  pressure  drop  relationships 


in  Chapter  III  and  a  control  strategy  in  Chapter  IV.  Finally,  the  results 
of  a  typical  computer  simulation  are  presented  and  discussed  in  Chapter  V. 
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II  Heat  Transfer  Analysis  of  Packed  Beds 

In  this  chapter,  the  study  of  the  flowrate  of  air  through  rock  beds  is 
examined.  Heat  transfer  coefficient  formulas  are  discussed  and  the  differ¬ 
ential  equations  which  describe  the  response  of  the  bed  are  developed.  A 
finite-difference  numerical  solution  to  these  equations  is  also  presented. 

Rock  Bed  Geometry 

The  rock  bed  is  a  relatively  simply  looking  storage  unit,  for  its  geo¬ 
metry  refer  to  Fig  1.  The  bed  consists  of  a  container  to  hold  the  rocks 
plus  an  inlet  and  outlet  for  the  flow  of  air  through  the  unit.  For  computer 
purposes  it  is  divided  into  N-equal  segments,  all  assumed  isothermal  and  of 
length  Ax;  this  assumption  becomes  better  as  the  number  of  segments  is  increased. 
Also,  insulation  surrounds  all  sides  of  the  bed  to  minimum  energy  losses. 

A  well  designed  bed  is  cubical  in  shape  and  utilizes  rocks  that  may  vary 
between  .5  and  2.5  inches  in  diameter;  these  two  facts  result  from  pressure 
drop  and  heat  transfer  considerations  (Ref  17). 

Rock  Bed  .Assumptions 

a.  Rocks  are  spherical  with  void  fractions  of  30  to  60  per  cent. 

b.  Fluid  must  be  a  gas  such  as  air. 

c.  Constant  fluid  and  material  properties  (for  air  this  is  a  very  good 
approximation  for  the  temperature  range,  70-200  F) . 

d.  For  practical  purposes ,  it  is  advisable  to  deal  with  a  uniform  or 
average  volumetric  heat  transfer  coefficient,  rather  than  focus  on  local  values. 

c.  Constant  volumetric  flowrates. 

f.  Fluid  flow  and  heat  conduction  in  the  solid  material  arc  one -d i mens ional 
in  the  flow  direction. 
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Fig  1.  Rock  Bed  Geometry 


g.  Thermal  storage  within  the  fluid  (air)  is  neglected. 

h.  Heat  conduction  within  the  fluid  (air)  is  neglected. 

i.  Fluid  properties  are  evaluated  at  a  film  temperature,  tp  where 

t  +  t 
_  o  “ 
tf - 5 - 


and 


t  =  free  stream  air  temperature 


t  =  rock  temperature 


Heat  Transfer  Coefficients 

There  are  several  methods  for  determining  the  heat  transfer  coefficient 
between  the  air  and  the  solid.  However,  in  this  section  a  relationship  which 
takes  into  account  both  laminar  and  turbulent  flow  contributions  will  be  con¬ 
sidered.  The  following  expressions  and  equations  were  cited  from  Ref  18. 

For  alternate  methods  refer  to  Appendix  B. 

The  void  fraction,  c  ,  is  defined  as 

g 


void  volume  of  bed  _  Vo 
c  g  total  volume  of  bed  V 


(4) 


therefore,  the  particle  or  rock  fraction,  e  ,  is 

_ rock  volume  \  - 

e n  total  volume  of  bed  1  <> 


5) 


For  a  bed  with  N  particles  which  have  a  volume,  V^,  the  volume  occupied  by 
the  rocks  is 


V  N  =  V 
P 


V(1  -  t  (J 


Then,  the  number  of  particles  per  imit  volimie  is 


N 

V 


(b) 


(1 


The  hydraulic  radius,  ,  is  defined  as 


P  _  void  volume  of  bed 
h  surface  area  of  particles 


=  _i£V 
A  X 
P 


(«) 


where  A  is  the  surface  area  of  a  particle.  Substitutin'’  for  ^  from  equation  (?) 
yields 


,,  _  Vp  s 

n  -iru-TjY 


9) 


For  heat  transfer  purposes,  a  characteristic  length  and  velocity  are  needed. 
Ref  IS  used  six  times  the  hydraulic  radius  for  the  characteristic  length,  !.*. 
Substituting  for  R^  from  equation  (9)  yields 


L*  =  6R,  = 


bVp  g 


A  1-ffi 

p 


10) 


A  spherical  particle  diameter,  I)  ,  is  defined  as 


D  =  (AW  A 
P  1  P 


01) 


where 


vp  =  up  'Vf> 
p 

A  =  nil  2 
P  P 


So  now  the  characteristic  length  becomes 


*  =  !>  r>- 
P  1-  A 


(12) 


The  characteristic  velocity,  u*,  is  defined  as  the  average  velocity  ol 
the  fluid  flowing  through  the  voids.  Prom  continuity 


8 


(13) 


pA  •  ,u* 
void 


jV  A 
o 


where  VQ  is  the  superficial  or  frce-strean  velocity.  The  void  area,  A^-^, 
is  related  to  A,  the  cross  sectional  area  of  the  bed  by 

Avoid  =  V  CH) 

Combining  co'iatf’on  (14)  on<j  equation  (13)  the  characteristic 
velocity  becomes 

✓ 

u*  =  VQ/Cg  (IB) 

This  is  also  referred  to  as  the  instititial  velocity. 

Finally,  an  expression  for  the  characteristic  Reynolds  number  can  be 

written  as 


Re*  =  u*  L*/v 


(16) 


where  v,  the  momentum  diffusivitv,  is 


v  =  ygc/p 


(17) 


Co^bin-'yvr  equations  (12),  (15)}(I7)  OO  yi  n 


Re*  =  DpG/ySc  (1  -  cg) 


(18) 


where  G  =  superficial  mass  velocity  =  PV  . 


The  characteristic  N’ussclt  number  is  defined  as 


Nu*  =  hb*/Kf 


(19) 


Substituting  for  L*  ft -an  ea:r;f  f''r.  (If)  Info  <  r. 


l 


Nu*  =  -JSfc 
Kf 


(20) 


where  h  is  the  average  coefficient  of  heat  transfer  and  Kp  the  thermal  con¬ 
ductivity  of  the  fluid. 

The  empirical  formulation  of  the  Nusselt  number  is  represented  by  the 
following  equation  (Ref  18)  : 

Nu*  =  (1.5  Re*1/2  +0.2  Re*2/3)  Pr1/3  (21) 


Reynolds  number  to  the  one-half  power  accounts  for  laminar  flow  effects  while 
7/3 

Re*  is  its  turbulent  counterpart . 

Combining  equations  (20)  and  (21)  and  solving  for  h  gives 


li 


.  u  -  b 

eg  Dp 


(0.5Re*1/2  +  0.2  Re*2//'3)  Pr1/3 


(22) 


where  Pr  is  the  Prandtl  number  of  the  fluid.  Remember,  all  fluid  properties 
arc  evaluated  at  a  film  temperature  as  discussed  earlier. 

liquation  (22)  is  valid  for  void  fractions  less  than  .65  and  Reynolds 
numbers  greater  than  or  equal  to  100.  This  equation  can  be  used  for  either 
spherical  or  cylindrical  shaped  objects,  and  it  can  be  modified  for  staggered 
tube  arrangements,  for  quick  approximations  of  h,  an  experimental  curve  of 
Nusselt  number  versus  Reynolds  number  is  provided  in  Ref  18. 

To  determine  the  volumetric  heat  transfer  coefficient,  hv,  the  following 
expression  is  used 

hv  =  hdA^/Adx  (23) 


whe  re 


A  =  cross  sectional  area  of  lied  segment 
Adx  =  total  volume  of  lied  segment 


10 


dA,,  =  total  surface  area  of  particles. 

Therefore, 

"> 

dAh  =  4ttRp“  [number  of  particles] 

=  4ttRp“  [Adx  (1  -  fg)/-Dp3/6] 
or 

d\  =  Adx  (1  -  €g)  (24) 

P 

After  substituting  equation  (?4)  into  equation  (22)  ,bv,bc  nones 

=  6h  (1  -  6g)/Dp  (25) 

Internal  Resistance 

To  account  for  internal  temperature  gradients  within  the  rocks  a  non- 
dimensional  heat  transfer  parameter  known  as  the  Biot  number  must  first  be 
described. 

The  Biot  number,  Bi,  is  defined  as 


Fig  2.  Spherical  Rock  Resistances 
where,  the  internal  resistance  is 


Ri  -  Rs/KsAs 


(27) 


and 

Ks  =  thermal  conductivity  of  solid 
As  =  surface  area  of  rock. 

Ro,  the  external  resistance  is 


Ro 


1 

7iAs 


(28) 


where  h  is  determined  by  equation  (22).  Combining  conations  (TC), 

(27)  and  (28)  yields 

81 '  -m- "  -®1-  (2!» 


It  is  only  when  the  Biot  number  is  less  than  0.1  can  heat  conduction  in 
the  solid  material  be  neglected.  In  this  case,  the  thermal  resistance  of  the 
solid  particles  is  small  in  comparison  to  the  thermal  resistance  of  the  con¬ 
vective  film  between  the  air  and  solid.  Therefore,  for  Bi  <  .1  the  effect 
of  temperature  gradients  within  the  solid  particles  can  be  ignored. 

A  Biot  number  correction  proposed  by  "Babcock”  in  an  article  written  by 
Jefferson  (Ref  9)  is  a  factor  which  can  be  incorporated  into  the  heat  transfer 
coefficient  to  account  for  the  following  effects:  axial  (static)  fluid-solid- 
fluid  conduction,  fluid  phase  axial  dispersion,  fluid-solid  heat  transfer 
resistance,  and  most  importantly,  internal  particle  thermal  resistance.  The 
result  according  to  Babcock  in 


hveff 


hv 

(1  +  BT757 


(30) 


where 


8  -  -^S  — 
P  ncs  +  HCA 


and 


HCS  =  heat  capacity  of  solid  *  psc  (1  “ 
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HCA  =  heat  capacity  of  air  =  Pacpa^g 

Note  that  HCA  is  much  less  than  HCS  since  the  density  of  air  is  much  less  than 
the  rock  density,  therefore  (?  is  approximately  one;  so  for  air  the  effect  of  s 
on  equation  (30)  can  be  neglected. 

Also,  the  formula  proposed  by  "Babcock"  was  proven  experimentally  for 
Biot  numbers  up  to  about  4. 

Development  of  Differential  Equations 

To  develop  the  equations  which  describe  the  response  of  a  rock  bed,  two 
energy  balances  on  the  following  control  volume  must  be  made.  (A  similar  set 
of  equations  can  also  be  found  in  Ref  6) .  The  first  energy  balance  is  for  the 
fluid  (air)  equation  while  the  other  is  for  the  solid. 


i  +  1 


•  • 


-^maCpa  Ta  +  ~  Ta  dx  -  E 

Ik  •—  ^ 


out 


Fig  3.  Differential  Bed  Segment 
The  time  rate  of  change  of  energy  within  the  control  volume  is 

3T 

E  stored  =  ^gpacpaAdx  — gp-  (31) 
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whore,  egAdx  is  the  void  volume.  Applying  the  conservation  of  energy  law 


dA, 

Ein  =  Eout  +  E stored  h  Adx  (Ta 


Tb)Adx 


(321 


where 

hdAi 

-  =  h^  as  seen  from  equation  (23). 

•  •  • 

With  some  rearrangement  and  substituting  for  E^n>  F.Qut  and  Fst;orej» 
equation  (32)  simplifies  to 


3T  3T 

-p  V  c  dx  -  £g  p  c  dx  -r^-  =  h  (T  -  T,  )  dx  (35) 

^a  o  pa  3x  &  a  pa  3t  vv  a  b 


Now,  dividing  through  by  p  c  £g  dx  and  substituting  u*  for  V  /eg  yields 

a .  pa  o 


3Ta  +  *  3Ta 
__  +  u*  - 


-K 


3x  p  c  tTa  '  V 
a  pa  6 


(34) 


But  since  the  heat  capacity  of  the  bed  (p  dvol  c  J  is  much  greater  than 

s  ps 

the  heat  capacity  of  the  air  (p  dvol  c  ),  the  temperature  of  the  air  as  a 

a  pa 

function  of  time  may  be  neglected.  Therefore,  equation  (34)  becomes: 


JiA 

macpa 


(T 


V 


(35) 


where  rii 


paA% 


u" 


The  following  control  volume  will  again  be  used  to  '.epresent  the  bed 

QL 

I^V^Tv/VJ 


temperature  as  a  function  of  time. 


Wvjyvw 


Fig  4.  Bed  Control  Volume 
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Energy  equations  for  this  case  are 


E .  .  ,  ,  =  m  c  [T  -  -  T  .  ] 

into  bed  a  pa  a , 1  a ,  1  ♦  i 

'  out  of  bed  =  ^1  ’  w^ere  Q^»  t*ie  of  energy  lost  from  the  bed  is 

Ql  =  Uj  (Tb  -  TOB) 


K.  A. 
ins  ms 


where 


=  conductivity  of  the  insulation 

A.  =  surface  area  of  the  bed 
ins 

t.  =  thickness  of  insulation  surrounding  bed 
ins  - 

and,  TOB  =  air  temperature  outside  the  bed.  Also, 


Stored  "  (I  ’  CS)pscps  Adx  w 


where;  (1  -  cg)  Adx  is  the  solid  volume. 


Therefore,  again  applying  the  conservation  law  and  rearranging  gives 


(1  -  cg)p  c  Aax  — —  =  m  c  [T  .  -  T  .  ]  -  Q, 

v  K'Ms  ps  3t  a  pa  L  a,i  a,i  ♦!  'L 


1  1) 

Solving  for  yields 


9T,  m  c  „  Qi 

D  ^  iX  pci  _ -  rrri  -  T  1  _  ' J _ 

W  =  (1  -  ‘■gJPj.c  aax"  1  a,i  a,i  +1  (1  -  Eg')rsc~  Aax 


Solution  of  Equations 

The  two  equations  now  derived  can  be  approximated  in  a  finite-difference 
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t 


form.  In  this  section  an  implicit  solution  scheme  will  be  presented.  This 
scheme  provides  stability  and  thus  allows  for  larger  time  steps,  therefore, 
computer  costs  are  reduced  while  still  maintaining  reasonable  results.  Refer 
to  Appendix  C  if  an  explicit  or  a  modified- implicit  scheme  is  desired.  For 
the  equations  that  follow^the  subscript  "i"  refers  to  the  distance  increment 
and  superscript  "j"  refers  to  the  time  step.  As  mentioned  before,  the  bed 
is  broken  up  into  N- equal  segments.  Therefore,  the  equations  are  used  to 
solve  for  N-unknown  bed  and  N-unknown  fluid  temperatures. 

Our  partial  differential  equations  in  an  implicit  form,  as  functions  of 
new  time  and  position,  become 


3TA  j-1 

3x  J 


-hy  A 

A  c 
a  pa 


(41) 


and 


i+1  = 


Ac  „ 
a  pa 


n - 


A  ax 


-  QL[a  -  eg)pscps  A4x] 


(42) 


Finite-differencing  equation  (41)  yields 


T  .  -  T  • 

a,t»i  a,i 

A  X 


i+1  _  -h„  A  L  j+l 


-\»- 
m  c 
a  pa 


-  T, 


a, i+1  b,i+l 


h  A  Ax 

After  some  more  rearrangement  and  letting  C.  =  - 

a  pa 


rV+1  =  *  T  j+1  +  -X-t  j 

y i+i  1+ci  a>1  1  +  b» 


i+1 


(43) 


(44) 


Since  is  always  positive,  the  quantity  (1/1  +  C^)  will  never  become 
negative;  therefore,  this  particular  form  of  the  air  equation  will  remain 
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stable  for  any  Ax.  In  both  the  explicit  and  modified- implicit  solution 
schemes.  Appendix  C,  a  stability  criteria  is  required. 

Equation  (42)  in  a  fully- implicit  form  becomes : 


T  3+l  . 

b.i+1  b,i+l  _  A  c 


T  j+1  _  ^acpa  T  j+1 
a,i  Q*  a,i+l 

•  ql/q* 


where:  Q*  =  (1  -  Gg)p  c  Aax. 


Multiplying  through  by  At  and  letting  C 2  = 


ft  c 

a  pa  At 


t  j+l  _'rj  +  C  |  T  i  + 

o.i+1  "  p,i+l  l2  I  a,i 


j+1  _  T  j+1 


,i+i  '  13*- 


This  form  of  the  bed  equation  will  always  be  stable,  for  any  At.  Mow 
substituting  for  from  equation  (46)  into  equation  (44)  and  doing  some 

more  manipulation  yields: 

T  M  =  1+C1C2  T  j+1 
a,i+l  l+Cj+C.jC2  a,i 


r2  xb,i+l  ‘  1+C1+C1C2  ~Q* 


where  =  Uj  (Ib^i+l  *  TOB)  as  in  equation  (37).  Equation  (47)  represents 
the  fully- implicit  finite-differenced  air  equation  as  a  function  of  time 
and  position  in  the  bed. 

Substituting  for  T from  equation  (47)  into  equation  (46)  yields 
the  final  bed  equation. 


tH  - 
lb,i+l  '  1+C 


_ T  j+1  +  ^C1  |"T  3 

1+C1C2  a,i  1+C-^q  b,: 


-  Qta. 
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Equations  (47)  and  (48)  will  completely  describe  the  response  of  the  rock  bed 
for  nodes  2  through  k. 

However,  in  the  discharging  mode  where  a  complete  reversal  of  the  bed 
temperature  distribution  is  required,  the  bed  temperature  at  node  1  must  also 
be  determined.  To  do  this,  a  second-order,  polynomial  curve  fit  is  used.  The 
curve  is  as  follows: 


+  BAX  +  CAx2  =  T,  V+1 
b2 ) 

+  2BAx  +  4CAx2  =  Tb^+1 
+  3BAx  +  9CAx2  =  Tb^+1 


(49) 


These  equations  were  solved  using  the  Gauss-Siedel  Reduction  method;  the 
resultant  bed  temperature  at  node  1  is 


(50) 


Since  monotonic  temperature  curves  were  usually  present  in  both  modes, 
equation  (50)  proved  to  be  a  fairly  accurate  extrapolation. 


Analytical  Solution 

An  analytical  result  serves  the  useful  purpose  of  evaluating  the  com¬ 
pleteness  and  validity  of  the  numerical  method. 

The  analytical  solution  proposed  by  Clark  and  Arpaci  in  Ref  5  is  as 
follows:  Tor  the  air  temperature 


VS’6)  ‘  Ta,o  =  D(0)F1  (S’6)  + 

D(4)  -  D(o) 

- 5 - Gj(S,«)  (51) 
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and  for  the  bed  temperature 


Tb(S,f)  "  Tb>0  =  D(o) f^  (S,5)  + 

D(«J  -  D(o) 

- - -  gi  (S,6 j  (52) 

51  1 


The  initial  and  boundary  conditions  are: 

T  (x,o)  =  T  =  initial  air  temperature 

a  cl  ,  O 

T^Cx^)  =  q  =  initial  bed  temperature 
T  (o,t)  =  T  (t)  =  inlet  air  temperature 


For  stability  purposes  these  equations  were  only  valid  for  small  time  steps 
and  small  distance  increments.  Also,  the  equations  were  derived  assuming 
negligible  heat  loss  from  the  bed.  In  equations  (51)  and  (52) 


DM  •  Ta(0)  -  Ta  o 


Wl>  ‘  W  -  Ta,o 


and 


S  = 


K  Pw 


Pc  A  .  , 
a  pa  void 


P  _  K  Pw  , 

5  '  p  c  A1 
s  ps 


('-#) 


(53) 

(54) 


where 


u*  =  interstitial  velocity 

\oid  =  flow  area 

A*  =  (1  -  £g)  A  =  rock  area 

t  ■  time  in  bed 
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x  =  position  in  bed 
h  =  average  heat  transfer  coefficient 
Pw  =  wetted  perimeter  of  rock 


and  Pw/A  =  wetted  perimeter;  but  since  R.  ,  the  hydraulic  radius,  can 
flow  area  n 

also  be  written  as  the  ratio  of  flow  area  to  wetted  perimeterjOne  obtains 
from  equation  (10) 


Pw/\oid  "  6/L*  <55> 

Subst ituting  for  L*  from  equation  (12)  into  equation  (55)  and  letting  A*  = 
dA^/Adx  from  equation  (24)  yields 

fV  Avoid  =  A*/Cg  (56) 


Also  note  that  in  Figs  5  and  6 
V  =  u* 


9  =  t 
1 


C 


1 


p  ps 

p  =  p 

•  a 


C  =  c 
P  pa 

A  =  A  .  , 
void 

The  functions,  ,  Gj  and  f^  and  g,  from  equations  (51)  and  (52)  can  be 
determined  graphically  as  functions  of  S  and  6  ,  or  they  can  be  obtained 
mathematically  with  the  use  of  Bessel  functions.  To  check  the  validity  of 
the  numerical  program,  a  constant  inlet  air  temperature  was  tested,  and  the 
graphical  method  was  employed. 


For  a  constant  inlet  air  temperature  case 

Ta(6n}  =  Ta(o);  n  =  2>  3* ■ * 

and 

D(6n)  =  D(o) 

Equation  (SI)  becomes 


'  X,  -  (T 

dLfi  ^  § 


(O)  - 


T«,,>F1 


(57) 


Equation  (57)  represents  the  analytical  air  equation  for  constant -inlet  air 
temperatures.  Equation  (52)  yields  the  analytical  bed  equation. 

Tb  -  V  fi  lV0)  •  T»..> 

For  the  graphs  of  F^  and  refer  to  Figs  5  and  6. 

An  analytical  result  (Ref  7)  for  the  bed  temperature  at  node  1  (x=o) 
can  be  obtained  from  the  explicit  bed  equation  in  Appendix  D,  namely: 


-hA*  <Tb  -  V 


integrating  and  rearranging  we  get: 


(58) 


; 

*b,o 


-hA*  } 

U  -  ^  l 


(59) 


At  x=o,  T  is  constant  for  a  particular  Ax,  therefore  equation  (59)  becomes: 

d 


mOb  -  t. 


L,0^] 


b  -  -hA* 

V  (1-  CgJPsCps 


At 


(60) 


and  again  rearranging  yields  : 


a,o 

a,o 


exp 


( 


-hA* 


T ~ 


At 


Scp^ 


(61) 
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where : 


Ty  o  =  initial  air  temperature 
Q  =  initial  bed  temperature 

liquation  (61)  represents  the  analytical  response  of  the  bed  at  x=o.  This 
result  can  be  used  for  comparison  purposes  with  the  numerical  curve  fit 
solution  at  node  1. 


^2 


,  dimensionless  gas  temperature 


( ^  =  1  signifies  no  heat  loss) 

Fig  5.  The  function  versus  <5  ^p*'  1  ^ 
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f,(S,8,ij}  DIMENSIONLESS  WALL  TEMPERATURE 


f 

This  chapter  deals  with  the  discussion  of  pressure  drop  relationships 
for  rock  beds,  flat  plate  collectors, and  duct  work. 


Power  Calculation 

Once  pressure  drop  has  been  determined  a  power  calculation  is  necessary 
to  obtain  fan  blowing  costs.  The  theoretical  work  per  unit  mass,  w^,  is 
given  by  the  pressure  drop  divided  by  the  air  density. 


w. 


AP 

Pi 


(62) 


Power,  the  rate  of  doing  work,  is  therefore 


P  =  power  =  wk*a  =  wk  (paAV) 


(63) 


In  terms  of  the  volumetric  flowrate,  Q  «  VA 


P 


(64) 


Power  estimates  can  be  easily  converted  to  operating  costs  by  using  standard 
utility  prices  and  a  fan  efficiency  factor. 


Rock  Bed 

Ergun's  equation  (Ref  3)  was  used  to  determine  rock  bed  pressure  drop. 
His  particular  relationship  is  valid  for  both  laminar  and  turbulent  flow's 
and  for  void  fractions  between  .40  and  .65. 


-AP 

L 


1-E 

170^ — S.  +  1.75 


(65) 


where 


p  =  density  of  fluid  (air) 
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V  =  free-stream  velocity 
o 

R  =  V  D  /vf  =  Reynolds  number 
ep  opt 

Ergun's  equation  can  also  be  written  in  the  following  form  which  accents  the 
presence  of  both  a  viscous  term  and  a  dynamic  pressure  term. 

-AP  =  apg  V  +  bpV2^  (66) 

°c  o  o 

where 

aug  V  =  viscous  term 
ac  o 

2 

bpV  =  pressure  term 


After  some  algebraic  manipulation 


A  tube  bundle  pressure  drop  expression  is  included  in  Appendix  B. 

Collector 

A  flat-plate  collector  consists  of  an  absorber  plate,  made  of  metal  for 
strength  and  good  thermal  conduction,  insulation  to  reduce  energy  losses 
from  both  the  bottom  and  sides  of  the  collector,  and  one  or  more  cover  plates 
of  either  glass  or  plastic  to  reduce  heat  losses  from  the  top  of  the  col¬ 
lector  (Fig  7) . 


t 


Collector  cross-section 


Cl 

ass  or  Plastic  Covori 

3 

1 

mbhkbhhmI 

in 

nwifii 

*  L»  Insulation 

Absorber  Plate 


Fig  7.  Collector  Geometry 

For  this  simple  flat-plate  air  heater: 

A  ^  ■  aB  =  air  duct  cross  sectional  area 
Pc  =  2  (a+b)  =  air  duct  perimeter 
a  =  B/a  =  aspect  ratio  of  collector 
L  =  length  of  collector 
B  =  width  of  collector 
a  =  collector  air  gap 

The  collector  pressure  drop  can  be  expressed  in  terms  of  a  friction 
factor,  f,  (Ref  10). 


A_  _  4fL pV2 


(69) 


where: 

f  =  Fanning-friction  factor 
p  =  air  density 

=  hydraulic-diameter 
Fbr  parallel  flat  plates: 

n  4AW  _  2aB 
°h  ~  P  a+F  ’ 


(70) 
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Substituting  for  in  equation  (69J  gives 

1P  ■  fit  <a*BJ  <7» 

*C 

For  fully  developed  turbulent  flow  in  tubes  (this  is  also  valid  for 
parallel  plates  if  the  hydraulic  diameter  is  used)  and  a  Reynolds  number 
range  of  5,000  to  30,000,  the  following  Fanning  friction  factor  may  be  used 
(Ref  10) : 

f  =  .079Re'-25  (72) 

and  for  30,000  <  Re  <  1,000,000: 

f  =  . 046Re" * 20  (73) 

where  Re  =  pVD^/yg^,  and  V  =  Q/Acol-  Now  substituting  equations  (72)  and 
(73)  into  equation  (71)  yields: 

(«) 

for  5,000  <  Re  <  30,000;  or 

4p  =  !p®‘’  v2(I^)'20(a*b,  t75) 


for  30,000  <  Re  <  1,000,000. 

The  Reynolds  number  for  collectors  is  normally  less  than  30,000  since  "a" 
is  usually  quite  small;  therefore  equation  (74)  is  frequently  used  in  calcu¬ 
lations.  This  equation  can  also  be  written  as 


AP 


.  079L 
gcP(aB)3 


(76) 
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A  D 

And  since  power  =  P  =  iti  — 


.25 


n  _  . 079L(a+b) f m  V  1 
®  p*  aB J  Re 


(77) 


Pressure  drop  is  extremely  sensitive  to  variations  in  "a''  or  the  collector 
air  gap  spacing.  As  "a"  increases,  AP  decreases;  but  since  ''a"  is  much 
less  than  "B",  the  increase  in  surface  area  for  heat  loss  is  minimal.  There¬ 
fore,  since  the  shape  of  the  air  duct  is  normally  arbitrary,  a  somewhat 
larger  gap  spacing  is  recommended  when  one  considers  the  significant  reduc¬ 
tion  in  pressure  drop  and  blowing  power  requirements. 

Duct  Work 

Our  design  analysis  consists  of: 

a.  straight,  round  smooth  radius  ducts 

b.  90°  -  elbows,  smooth  and  round. 

For  most  solar  energy  applications,  the  duct  design  is  fairly  simple,  involving 
mainly  straight  channels  and  elbows,  therefore,  pressure  drop  calculations 
will  be  considered  for  these  two  important  and  essential  cases  only. 

90°  elbows:  To  determine  pressure  drop  in  elbows  one  needs  to  know  the 
loss  coefficient  which  is  represented  by  the  ratio  of  the  total  pressure  loss 
to  dynamic  pressure  at  a  referenced  cross  section,  "oM  usually  at  the  beginning 
of  the  configuration. 


or 


C  =  loss  coefficient 
o 


(78) 


(79) 


ye 


X 


J 


where  "C^"  can  be  found  in  the  ASHRAE  Handbook  of  Fundamentals.  The  geometry 
for  elbows  as  depicted  in  ASHRAE  (Ref  1)  is  shown  in  Fig  8. 


Fig  8.  Duct  Elbow 


C  =  KC  1  (80) 

o  o  1 

where  C  *  is  dependent  on  r/D  ratios  and  K  is  dependent  on  e. 

Note  that  K  *  1  when  e  =  90°. 


Coefficients 

for  90°  Elbows 

r/D 

0.5 

0.75 

1.0 

1.5  2.0 

2.5 

CoX 

0.71 

0.33 

0.22 

0.15  0.13 

0.12 

For  angles  other  than  90°  multiply  by  the  following  factors: 


0  0  20  30  45  60  75  90  110  130  150  180 

K  0  0.31  0.45  0.60  0.78  0.90  1.00  1.13  1.20  1.28  1.40 


For  design  purposes  each  straight  channel  is  referred  to  as  a  "length" 


and  each  90°  elbow,  a  "turn".  The  motivation  for  circular  ducts  as  opposed 
to  rectangul ar  ducts  is : 


a.  simpler  to  manufacture 

h.  less  leakage  and  easier  to  insulate 

c.  less  costly 

d.  most  important;  less  surface  area,  therefore,  less  pressure  drop 
and  thus  less  blowing  power  required. 

Therefore,  for  these  reasons,  it  should  prove  likely  that  ones  final 
decision  would  be  to  use  circular  cross  sections.  The  only  real  advantage 
with  rectangular  ducts  is  probably  that  they  fit  neatly  in  between  beams, 
thus  creating  somewhat  of  a  hideaway  design. 

Straight  lengths:  Pressure  drop  in  straight  channels  can  again  lie 
expressed  in  terms  of  the  Fanning  friction  factor,  f. 


-AP  = 


4  fhp\'“ 


181) 


where  Dj  =  - =  Ilj, ;  or  simply  the  actual  tube  diameter  when  dealing  with 
circular  ducts. 

As  before, 

()7q 

f  =  for  5,000  <  Re  <  30,000 

Re 

and  f  =  for  30,000  <  Re  <  1,000,000. 

Re 


Therefore, 


(82) 


(S3) 
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where : 


L  =  duct  length 
Re  =  pVDh/ugc 

After  substituting  equation  (82)  into  equation  (63) 


(84) 


A  rule  of  thumb  (ASHRAIi,  Ref  1)  for  duct  work  is  that  duct  velocity  should  not 
exceed  700  fpm  for  low  noise  purposes; 


Duct  velocity  =  V  = 


volumetric  flowrate  _  Q 

cross  sectional  area  of  ducf  -  A, 


D 


where:  A^  =  ^1)^,  / 4 


Design  recommendations  for  pressure  drop  in  all  three  configurations  arc 
(Ref  2) : 

a.  rock  beds  -  .1  to  .4  in  H?o 

b.  collectors  -  .2  to  .8  in  ll0o 

c.  duct  work  -  .08  in  H^o/100  ft  duct  length 

Note  that  duct  heat  losses  had  to  lie  minimal  to  validate  the  assumption 
that  the  inlet  air  temperature  for  the  bed  during  its  charging  mode  was  nearly 
the  same  as  the  air  temperature  out  of  the  collector.  This  involved  the 
prediction  of  heat  transfer  from  ducts  as  presented  and  discussed  in  Appendix  ( 


IV  Performance  Simulation  of  an  Air  Heating  System 

This  chapter  is  concerned  with  developing  a  model  to  predict  the  oper¬ 
ational  modes  of  an  air  heating  system.  Also,  a  space  heating  load  model 
along  with  auxiliary  energy  requirements  are  discussed. 

Control  Strategy 

A  suitable  method  for  long  term  simulations,  is  to  assume  that  during 
any  time  period,  the  system  operates  in  whatever  modes  necessary  to  maintain 
the  building  temperature  at  the  desired  level.  During  each  time  period,  the 
amount  of  energy  collected  is  compared  with  the  amount  of  energy  required  to 
meet  the  load  (Ref  13) . 

To  fully  complete  a  simulation  model  there  are  three  modes  that  must  he 
satisfied  for  a  solar  air  heating  system.  For  the  workings  of  an  air  system 
refer  to  the  schematic  diagrams  in  Figs  9  and  10,  note  that  only  one  fan  is 
necessary  for  its  operation. 

In  the  first  mode  of  operation  solar  energy  is  available  but  not  enough 
is  collected  to  meet  the  space  heating  load;  therefore,  the  collected  energy 
is  delivered  directly  from  the  collector  to  the  house.  In  mode  two,  solar 
energy  is  sufficient  to  meet  the  space  load;  therefore,  the  energy  collected 
is  delivered  to  the  house  until  the  load  is  satisfied;  then  for  the  remaining 
time  period  the  rest  of  the  energy  is  transferred  to  the  storage  bed.  This 
is  commonly  referred  to  as  the  charging  mode,  flic  third  mode  becomes  operational 
when  solar  energy  is  no  longer  available  (usually  at  night);  then  hot  air  is 
drawn  from  the  hot  end  of  the  bed  until  the  space  load  is  met  and  room  temperature 
air  is  returned.  This  is  known  as  the  discharging  mode.  In  modes  one  and 
three,  auxiliary  energy  may  be  required  if  the  amount  of  energy  transported 
from  cither  the  collector  or  the  storage  unit  is  depleted  before  satisfying 
the  load. 
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Schematic  Diagram  of  a  Solar  Air  Heating  System  (Ref 
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Fig  10.  Three  Modal  Diagram  of  Control  Box 


A  variety  of  factors  influence  space  heating  requirements  such  as  the 
location  of  the  building,  its  design,  orientation  and  particular  habits  of 
the  occupants. 

Many  space  heating  models  have  been  proposed,  ranging  in  detail  from 
simple  to  relatively  complex  approaches.  However,  due  to  computer  time 
and  costs,  complex  models  for  load  calculations  are  not  justified.  There¬ 
fore,  a  simple  space  heating  load  model,  the  degree-day  model,  expresses 
the  heating  load,  QSHL,  as  a  function  of  the  difference  between  the  inside 
temperature  of  the  building  (reduced  to  account  for  heat  generation  due  to 
lights,  people,  etc;  usually  about  651')  and  the  mean  ambient  (outside)  tem¬ 
perature  (Ref  12). 

qm  -  UA  (Tinsidc  -  Ta)  (85) 

where : 

UA  =  the  building  overall  loss  coefficient-area  product 
The  number  of  degree-days  in  a  month,  DD,  is 

N 

DD  =  T.  (651-  -  T  )  (So) 

i=l  a  1 

where  N  is  the  number  of  days  in  a  month  and  (6SF  -  T  ).  is  taken  to  be  _ero 

a  i 

for  mean  daily  temperatures  above  65F. 

Therefore,  a  monthly  space  heating  load,  QSIIIA1,  can  be  written  as 

QSIILM  =  UADI)  (8") 


For  simulation  purjxjscs,  an  hourly  heating  load  was  used,  namely: 


Auxiliary  energy  is  provided  for  space  heating  whenever  the  amount  of 
solar  energy  collected  is  insufficient  to  meet  the  space  heating  load.  In 
practice,  this  condition  is  detected  by  a  thermostat  monitoring  the  temper¬ 
ature  within  the  building.  In  the  present  space  heating  model  however, 
the  temperature  within  the  building  is  assumed  to  be  constant  and  auxiliary 
energy  is  supplied  whenever  the  rate  at  which  solar  energy  can  be  provided, 
QAVG,  is  less  than  the  space  heating  lond,QR,rT,,  "’herefore,  the  amount  of 
auxiliary  energy,  nATTX,  needed  is: 

PATTYsOgpr  _  oyprp  _  OT.oom 

where  it  has  been  assumed  that  any  hed  ’’osses  (ot  eo-p  actua"!‘,v  re-him 
space  heating  load  (Ref  IP), 

In  the  discharging  mode,  when  energy  Is  tairer.  4’v,om  4_Ke  storage  unit , 
the  auxiliary  energy  needed  is: 

QA’JX  =  Qf!rr  —  OPT7T''H  (op) 

where, QR,rnR,  the  rmiount  of  energy  drawn  from  the  bed  is: 

QRRIHH  =  mcpa(TFB  -  'T’ROTiM)  (01) 

and 

rnr,".R  =  Rxit  air  temporal  lire  of  bed 

TROOM  =  Room  air  temperature 

^'b.o  QATTV.  calculations  were  used  to  determine  backup  energy  costs  In  the 
air  system  simulation. 

^or  the  ccmnlet e  logic  of  an  air  system  Simula4  Ion,  a  owchor4  e4'  tpe 
computer  program  Is  presented  in  ''nnendlx  77 . 


V  Results  and  Conclusions 


♦ 

9 

Rock  Bed  Simulation  Results 

The  major  objective  in  using  a  fully-implicit  finite-difference  solution 
scheme  was  to  be  able  to  use  time  steps  as  large  as  one  hour  and  still  obtain 
adequate  results  for  realistic  simulations.  In  order  to  justify  this  criteria 
the  rock  bed  was  tested  separately  for  a  constant  inlet  air  temperature  and 
specified  initial  conditions  as  indicated  in  Table  1. 

For  this  particular  rock  bed  simulation,  the  air  and  bed  temperatures 
are  somewhat  more  sensitive  to  changes jn Ax than  to  changes  in  time  step.  The 
corresponding  temperatures  for  a  Ax  of  six  inches  and  a  At  of  one  hour  were 
reasonably  close  for  engineering  purposes  to  the  temperatures  for  a  one  inch 
distance  increment  and  a  .1  hour  time  step.  Therefore,  these  particular 
increments  were  used  in  the  air  system  simulation 

To  further  appreciate  this  selection,  the  air  and  bed  temperatures  were 
plotted  as  functions  of  position  and  time  in  the  bed  at  various  time  steps. 

In  Fig  11  one  notices  that  for  a  At  of  one  hour  the  air  temperatures  near  the 
beginning  of  the  bed  were  slightly  lower  than  those  corresponding  to  a  At  of 
.1  hour.  However,  at  x=2  feet  the  trend  began  to  reverse  itself  as  it  should, 
since  energy  must  be  conserved.  But  for  all  positions  in  the  bed,  the  temper¬ 
ature  difference  between  the  largest  and  smallest  time  increments  were  again 
relatively  small.  In  Tig  12,  the  same  general  trend  is  observed.  Only  for 
this  case,  the  bed  temperatures  arc  a  bit  lower  than  the  corresponding  air 
temperatures  of  Fig  11,  but  as  the  volumetric  heat  transfer  coefficient  approaches 
infinity,  the  air  and  bed  temperatures  at  a  particular  point  in  t he  bed  become 
nearly  equal.  Note  that  in  both  these  figures  a  Ax  of  1  foot  was  used;  therefore, 
a  smaller  Ax  would  obviously  result  in  reducing  the  temperature  difference 
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Table  1.  Temperature  Convergence  Data 

Case  Tested:  Rock  Diameter  =  1  inch 

Initial  Bed  Temperature  =  70F 
Inlet  Air  Temperature  =  140F 
Volumetric  Floivrate  =  1200  CFM 

B  2 

Average  Heat  Transfer  Coefficient  =  3.22  /hr-ft“-F 


Ax  (inches) 

At  (hours) 

*TA  (F) 

*TB  (F) 

1.0 

0.1 

96.1 

91.1 

1.0 

0.5 

94.94 

91 . 3 

1.0 

1.0 

93.37 

90.5 

4.0 

1.0 

96.5 

93. 3 

**6.0 

1.0 

98.2 

94.7 

*  -  air  5  bed  temperatures  at  x=l  foot  after  one  hour  of  charging 
**  -  increments  used  for  the  simulation  runs 
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Temperature  Profile 


between  the  va rious  time  steps.  This  again  reinforces  the  choice  of  selecting 
a  Ax  of  six  inches  and  a  time  step  increment  of  one  hour. 

In  Fig  13,  the  air  temperatures  were  plotted  as  a  function  of  time  for  a 
particular  x- location  in  the  bed  (x=l  foot).  With  increasing  time,  the  air 
temperatures  at  x=l  foot  began  approaching  140F,  the  inlet  air  temperature,  TIN. 
Also,  the  air  temperatures  for  the  one  hour  time  step  were  all  slightly  lower 
than  the  corresponding  air  temperatures  at  the  other  three  time  increments. 

These  two  occurrences  followed  both  a  natural  and  expected  trend. 

To  test  the  validity  of  the  numerical  scheme  an  analytical  solution  as 
discussed  in  Chapter  II  was  used  for  comparison  purposes.  The  analytical 
solution  was  limited  to  relatively  small  time  and  distance  increments,  for 
stability  purposes.  The  solution  also  assumed  negligible  heat  loss  from  the 
bed;  this  was  simulated  in  the  numerical  case  by  simply  setting  the  conductivity 
of  the  insulation,  KINS,  to  zero.  In  Figs  14  and  15,  the  air  and  bed  temper¬ 
atures  were  plotted  as  functions  of  x  for  a  charging  period  of  only  .1  hour. 

In  both  graphs  the  numerical  solution,  represented  by  the  circles,  was  extremely 
close  to  the  analytical  result.  Also  note  that  for  the  first  few  inches  of  the 
bed  the  air  temperatures  were  significantly  larger  than  the  corresponding  bed 
temperatures,  and  a  x=12  inches,  both  the  temperatures  approached  70F,  the 
initial  bed  temperature. 

In  Figs  16  and  17,  the  air  and  bed  temperatures  were  plotted  as  functions 
of  time  at  x=l  inch,  for  a  At  of  .1  hour.  Again,  the  numerical  temperatures 
were  quite  similar  in  comparison  to  the  analytical  ones.  Also  for  increasing 
time  the  air  and  bed  temperatures  approached  140F,  the  constant- inlet  air 
temperature;  at  t=l  hour  the  air  temperature  was  approximately  139F  while  the 
corresponding  bed  temperature  was  around  138F. 
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In  addition,  in  Fig  15  the  numerical  curve  fit  extrapolation  at  x=o 
resulted  in  a  bed  temperature  of  98F;  and  the  dotted  line  joining  the  bed 
temperatures  at  x=l  and  x=o  obviously  displayed  a  rather  nice  consistency 
with  the  natural  extension  of  the  solid  analytical  curve. 

Therefore,  Figs  14  -  17  essentially  prove  the  validity  of  the  numerical 

scheme  and  also  the  second-order  polynomial  curve  fit  solution  at  x=o. 

The  next  three  figures,  Figs  18,  19  and  ZO^epresent  the  bed,  collector 
and  duct  work  pressure  drops,  respectively.  The  bed  pressure  drop  is  plotted 
as  a  function  of  rock  diameter  for  various  flowrates,  while  the  collector  and 
duct  work  pressure  drops  are  plotted  as  functions  of  Q,  the  volumetric  flowrate. 
It  can  be  seen  from  Fig  18  that  the  bed  pressure  drop  varies  directly  with  Q 
and  is  inversely  proportional  to  the  rock  diameter.  For  rock  diameters  less 
than  one  inch  significant  pressure  drop  may  result,  especially  at  the  higher 
flowrates.  In  Figs  19  and  20,  both  the  collector  and  the  duct  work  pressure 
drop  vary  directly  and  almost  linearly  with  the  flowrate.  Here,  it  is  worth¬ 
while  to  mention  that  even  though  heat  transfer  coefficients  increase  with  Q. , 
pressure  drop  and  thus  blowing  costs  also  increase.  However,  this  cost  increase 
generally  outweighs  any  possible  gain »» -tin thermal  performance  of  the  rock  bed 
that  results  due  to  increasing  flowrates.  Therefore,  one  is  normally  advised 
to  use  existing  rules  of  thumb  which  provide  the  satisfactory  range  of  flowrates 
per  square  foot  of  collector  area. 

The  last  two  figures  represent  total  monthly  blowing  costs  as  related  to 
various  rock  diameters.  Total  cost  refers  to  the  combination  of  bed,  collector 
and  duct  work  blowing  costs.  Monthly  costs  arc  also  directly  related  to  Q  but 
inversely  related  to  rock  diameter.  Blowing  costs  remain  relatively  constant 
for  rock  diameters  of  one  inch  or  more.  However,  for  larger  diameter  rocks 
the  rate  of  heat  transfer  is  substantially  reduced,  since  the  heat  transfer 
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coefficient  is  inversely  proportional  to  rock  size.  Therefore,  these  two 
considerations  alone  appear  to  suggest  that  rock  diameters  of  approximately 
one  inch  should  prove  most  economical.  According  to  Fig  21,  for  one  inch 
rocks,  the  blowing  costs  at  1200  CFM  correspond  to  approximately  two  dollars 
per  month  and  at  1800  CFM,  about  six  dollars  per  month;  these  results  were 
based  on  an  eight  hour  period  of  operation. 

In  the  air  system  simulation,  a  volumetric  flowrate  of  1200  CFM  with  a 
bed  of  one  inch  diameter  rocks  was  used.  Therefore,  for  this  particular 
simulation  one  could  safely  assume  that  the  total  operating  costs  of  two 
dollars  per  month  contributes  rather  insignificantly  when  determining  the 
life  cycle  cost  effectiveness  of  the  air  heating  system.  Also,  with  the 
above  conditions  but  for  24  hours  of  operation,  the  resulting  blowing  cost, 
as  indicated  in  Fig  22,  is  approximately  $5.50  per  month;  this  however  is  still 
a  relatively  small  amount  when  one  considers  the  total  investment  costs  required 
per  square  foot  of  collector  area. 
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BED  TEMP  VERSUS  BED  POSITION 


Comparison  of  Analytical  and  Numerical 
Bed  Temperatures  (function  of  x) 


Comparison  of  Analytical  and  Numerical 
Air  Temperatures  [function  of  time) 


Comparison  of  Analytical  and  Numerical 
Bed  Temperatures  (function  of  time) 
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Fig  18.  Bed  Pressure  Drop 


COLLECTOR  PRESSURE  DROP  VERSUS  FLOW  RRTE 


Collector  Pressure  Drop 
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System  Simulation  Results 

For  illustrative  purposes  the  computer  program  was  used  to  simulate  a 
residence  in  Columbus,  Ohio.  Data  for  the  weather,  position  and  orientation 
of  the  collector  is  the  same  as  that  used  in  Prins'  program  (Ref  16).  The 
parameter  values  used  are  as  follows: 

1.  Latitude  =  40  degrees 

2.  Number  of  covers  =  2 

3.  Slope  of  collector  =55  degrees 

4.  Azimuth  angle  =  0  degrees 

5.  Area  of  collector  (varied  parameter)  =  400,  500,  600  and  700  sq.  ft. 

6.  Thickness  of  cover  material  =  .23  cm 

7.  Extinction  coefficient  of  cover  material  =  .161/cm 

8.  Index  of  refraction  of  cover  material  =  1.526 

9.  Absorptivity  of  the  absorber  plate  =  .92 

10.  Heat  removal  efficiency  =  .62 

11.  Heat  loss  coefficient  =  .60 

12.  Mass  flowrate  of  fluid  through  collector  =  5112  lbm/hr 

13.  Specific  heat  of  collector  fluid  =  .24  B/lbm-F 
1',.  Dirt  factor  =  .02 

15.  Initial  temperature  of  air  in  bed  =  70. OF 

16.  Efficiency  of  heat  exchanger  =  1.0 

17.  System  in  use  =4.0 

18.  Number  of  months  program  is  to  run  =  8 

19.  Time  increment  for  discharge  =  2.0  minutes 

20.  Epsilon  for  time  increment  =  .01  hours 

21.  Minimum  pump  energy  =  275.0  B/hr 

22.  Structure  conductance  =  1000.0  B/hr-F 


55 


23.  lime  step  =  1.0  hours 

24.  Volumetric  flowrate  =  20.1  cfs 


25.  Bed  height  =  8.0  ft 

26.  Bed  width  =  10.0  ft 

27.  Air  density  =  .071  lbm/ft^ 

28.  Rock  diameter  =  .08333  ft 

29.  Prandtl  number  for  air  =  .711 

30.  Conductivity  of  air  =  .0155  B/hr-ft-F 

2 

31.  Momentum  diffusivity  of  air  =  .649  ft  /hr 

32.  Specific  heat  of  air  =  .24  B/lbm-F 

33.  Air  temperature  outside  bed  =  70. OF 

34.  Rock  density  =  165  lbm/ft^ 

35.  Specific  heat  of  rock  =  .21  B/lbm-F 

36.  Bed  length  =  10.0  ft 

37.  Number  of  nodes  in  bed  =  21 

38.  Void  fraction  of  bed  =  .42 

39.  Conductivity  of  rock  =  1.0  B/hr-ft-F 

40.  Initial  bed  temperature  =  70. OF 

41.  Conductivity  of  insulation  =  .023  B/hr-ft-F 

42.  Thickness  of  insulation  =  1.0  ft 

2 

43.  Surface  area  of  insulation  =  420.0  ft 

2 

44.  Dynamic  viscosity  of  air  =  .046  lbm-hr/ft 

45.  Fan  efficiency  =  .55 

46.  Collector  air  gap  spacing  .3333  ft 

47.  Width  of  collector  =  4.0  ft 

48.  Length  of  collector  =  8.0  ft 

49.  Length  of  duct  work  =  200.0  ft 


50.  Duct  diameter  =  1.5  ft 

51.  Number  of  elbows  =  5.0 

52.  Hours  of  operation  =  8.0  hours 

53.  Primary  Federal  tax  credit  =  .40 

54.  Secondary  Federal  tax  credit  =  .30 

55.  Ohio  tax  credit  =  .10 

56.  Initial  tax  credit  amount  =  2000.0  dollars 

57.  Secondary  tax  credit  amount  =  1000.0  dollars 

58.  Maximum  tax  credit  amount  =  4000.0  dollars 

59.  Annual  mortgage  interest  rate  =  .12 

60.  Down  payment  =  .10 

62.  Backup  furnace  cost  =  $3. 26/million  Btu  (Gas) 

63.  Conventional  furnace  cost  =  $3. 26/million  Btu  (Gas) 

64.  Efficiency  of  solar  backup  furnace  =  .55 

65.  Efficiency  of  conventional  furnace  =  .55 

66.  Property  tax  rate  =  0.0 

67.  Income  tax  bracket  =  .40 

68.  Extra  insurance  and  maintenance  costs  =  .01 

69.  General  inflation  rate  =  .08 

70.  Fuel  inflation  rate  =  .12 

71.  Discount  rate  (varied  parameter)  =  .10  and  .20 

72.  Term  of  economic  analysis  =  20.0  years 

73.  Area  dependent  costs  (air  collector)  =  $15/sq.  ft 

74.  Area  independent  costs  (fans,  etc.)  =  1000.0  dollars 
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Three  parameters  were  varied  in  this  simulation  while  the  others  were 
held  fixed.  The  varied  parameters  were  the  collector  size,  the  discount 
rate  and  the  mortgage  loan  term.  Varying  the  collector  area  was  necessary 
to  determine  the  optimum  collector  size  for  maximum  savings.  A  cost  per¬ 
formance  curve  is  shown  below. 


Net 

Savings 
($)  9000 


Optimum  collector  size 


Fig  23.  Net  Savings 

Net  savings  refers  to  the  total  fuel  savings  minus'  investment  costs. 


58 


For  this  simulation,  the  optimum  collector  area  was  approximately  600 
square  feet,  this  corresponded  to  a  maximum  savings  of  8572  dollars  for  a 


discount  rate  of  10  per  cent  and  a  20  year  mortgage  term.  At  a  discount 
rate  of  20  per  cent  a  slight  savings  was  still  achieved,  therefore  the 
system  provided  at  least  a  20  per  cent  return  on  investment.  In  this  case, 
it  is  interesting  to  note  that  the  amount  of  savings  actually  increases 
with  lon^ct-  mortgage  periods.  This  is  due  to  the  fact  that  the  rate  of 
return  (20%)  is  greater  than  the  interest  rate  for  the  mortgage  loan.  It 
is  also  important  to  note  that  fuel  costs  were  based  on  natural  gas  prices , 
the  cheapest  of  the  fossil  fuels. 

For  other  key  results  of  our  simulation  refer  to  Tables  2  and  7>.  A  plot 
of  solar  energy  percentage  versus  collector  area  is  presented  below.  Note 
that  the  slope  of  the  curve  begins  to  drop  off  at  higher  collector  areas. 
Because  a  substantial  portion  of  the  solar  system  cost  increases  linearly 
with  collector  area  and  thermal  performance  increases  less  than  linearly, 
a  maximum  amount  of  savings  results  at  some  optimum  collector  size. 
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Fig  24.  Per  Cent  Solar  F'ncrgv 
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Collector  area  (ft^) 
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Structure  Conductance  =  1000  B/hr- 
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General  Remarks  on  Rock  Red  Air  S\ steins 

The  basic  advantages  of  a  solar  air  heating  system  are  that  the  heat 
transfer  fluid,  air,  will  not  freeze  or  boil  and  materials  used  in  the 
system  generally  have  a  long  life  expectancy .  Maintenance  of  the  system 
is  expected  to  be  low  if  quality  collectors  are  used,  and  corrosion  problems 
are  usually  quite  minimal  when  compared  to  water  systems.  Normal  maintenance 
is  confined  to  cleaning  or  replacement  of  filters  and  proper  care  of  the 
blowers  or  fans.  Another  advantage  of  the  air  system  is  that  the  house  can 
be  heated  directly  from  the  collector,  the  mode  1  operation. 

The  rock  bed  storage  lias  the  following  advantages  (Ref  17): 

a.  Rocks  are  plentiful  and  cheap  to  purchase. 

b.  They  are  non-toxic,  non-corrosive  and  non-flammable. 

c.  Rocks  have  a  high  heat  storage  capacity. 

d.  Thermal  energy  is  efficiently  transmitted  to  the  rocks  by  air  cir¬ 
culating  through  the  bed  due  to  the  large  heat  transfer  area.  Also,  heat 
conduction  through  the  bed  of  rocks  is  low,  because  the  area  of  contact 
between  the  rocks  is  small;  this  leads  to  low  heat  losses  from  the  bed. 

These  advantages  make  the  rock  bed  a  fairly  efficient  storage  unit  because 
air  can  leave  the  bed  at  a  temperature  nearly  equal  to  that  of  the  rocks  at 
that  point,  which  in  turn  is  nearly  equal  to  the  temperature  of  the  hot  air 
entering  the  bed. 

The  disadvantages  of  a  rock  bed  air  system  arc: 

a.  Significant  blowing  power  may  result  at  high  volumetric  flowrates, 
for  water  systems,  pumping  power  requirements  arc  usual lv  minimal. 

b.  Rock  beds  generally  require  a  larger  volume  for  thermal  storage  than 
water  storage  tanks. 

For  design  purposes  a  rule  of  thumb  of  1  to  -1  clin  per  square  foot  of 
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collector  area  is  recommendei(Ref  2).  Also,  the  economic  optimum  storage 
capacity  is  0.50  to  2.00  cubic  feet  of  rock  storage  per  square  foot  of 
collector  are  (Ref  17).  Washed  stones  or  crushed  rock  one  half  inch  to  two 
inches  in  diameter  should  be  used.  Die  type  of  rock  does  not  affect  the 
performance  of  the  storage  unit  as  long  as  the  density  is  approximately 
160  lbm/t't'  (Ref  17).  There  are  two  kinds  of  rocks:  igneous  rocks  and 
sedimentary  rocks.  Table  4  shows  the  various  densities  of  these  two  types 
of  rocks.  Also,  Table  S  provides  a  list  of  some  common  heat  storage 
materials  (Ref  17). 

Conclus ions 

In  conclusion,  a  solar  air  heater  can  effectively  provide  space  heating 
in  residential  buildings,  f.ven  when  natural  gas  prices  were  used,  the  air 
system  still  remained  cost  effective.  Tor  this  particular  simulation,  greater 
than  20  per  cent  rate  of  return  on  investment  was  achieved.  The  numerical 
method  proved  to  be  stable  and  convergent  and  showed  satisfactory  agreement 
in  comparison  to  an  analytical  solution  for  constant -inlet  air  temperatures. 
.Also,  the  three-modal  simulation  model  demonstrated  the  capability  of  cstimatin 
the  long  term  response  of  both  the  air  system  and  the  rock  bed.  In  addition, 
blowing  costs  proved  to  be  relatively  insignificant  in  this  particular  simu¬ 
lation,  resulting  in  a  total  operating  cost  of  only  about  two  dollars  per 
month.  Lastly,  computer  costs  were  relatively  low,  approximately  one  dollar 
per  simulation. 
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Table 

4 .  Types 

of  Rocks  and  Their  Densities 

IGNEOUS  ROCK 

SEDIMENTARY  ROCKS 

Name 

Density 
(lbm/ft  j 

Name 

Density 
(lbm/ft‘ ) 

Granite 

163-172 

Limestone 

162-169 

Syenite 

163-172 

Dolomite 

162-169 

Ryolite 

163-172 

Sandstone 

156-169 

Trachyte 

160-172 

Conglomerate 

155-169 

Diorite 

175-181 

Quartz  diorite 

169-178 

Dacitc 

169-178 

Andesite 

169-178 

Basalt 

181-194 

Obsidian 

144-169 
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Table  5.  Properties  of  Some  Cornnon  Heat  Storage  Materials 


Substance 

Specific  Heat 
cp 

(B/lbm-F) 

Densitv 

P  5 
(lbm/ ft' ) 

Thermal  Diffusivitv 

Air 

.240 

.075 

.850 

Asbestos 

.250 

9.40 

.032 

Alcohol 

.55-. 65 

50.0 

.0035 

Brine 

.71 

75.0 

.0100 

Brick 

.200 

106.0 

.0128 

Charcoal 

.242 

25.0 

.0660 

Concrete 

.200 

137.0 

.320 

Glass 

.16-. 20 

157.0 

.174 

Gasoline 

.  550 

50.0 

.0031 

tee 

.500 

57.3 

.1132 

Porcelain 

.255 

144.0 

.0160 

Rock 

.210 

160.0 

.550 

Sand  (dry) 

.191 

100.0 

.0100 

Steel 

.120 

490.0 

.482 

Water 

1.00 

62.4 

.0054 

Wood  (oak) 

.570 

48.0 

.0044 
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Appendix  A 


Review  of  Method  for  Determining  Useful  Fnergv  Collected 

The  collector  heat  removal  factor,  Fp,  of  equation  (1)  from  Chapter  II 
is  defined  as  a  quantity  that  relates  the  actual  useful  energy  gain  of  a  col¬ 
lector  to  the  useful  gain  if  the  whole  collector  surface  were  at  the  fluid 
inlet  temperature.  Fp  is  represented  by  (Ref  2) 


mcc  r(Tr  -Tr  .) 
f  pP  f,o  fjj/ 


FR  [(rajT^  -IJ 


(Tf,i'^ 


a” 


CA1) 


where : 

(~)fTt  =  amount  of  solar  energy  absorbed  by  collector  plate 
U  =  energy  loss  coefficient  of  collector 
The  effect  of  FR  is  to  reduce  the  calculated  useful  energy  gain  from  what 
it  would  have  been  had  the  whole  collector  been  at  Tf  -  to  what  actually 
is  using  a  fluid  that  increases  in  temperature  as  it  flows  through  the 
collector  (Ref  2).  As  the  mass  flowrate  through  the  collector  increases, 
the  temperature  rise  through  the  collector  decreases. 

Tp  in  equation  (Al)  is  given  by 

ITt-RrtH  (A2) 


where  r  is  the  ratio  of  hourly  total  radiation  on  a  horizontal  surface  to 
daily  total  radiation  on  a  horizontal  surface.  A  much  better  approximation 
of  r  than  the  one  Prins  used  has  now  become  available  (Ref  14).  The  new 
equation  is*. 


r 


t 


tj  (a+b  cos  w) 


cos  w  -  cos  ws 
s i n  ws  -  ws  cos  ws 


(AS.) 
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where 


M 


a  =  .409  +  .5016  sin  (ws-1.047) 
b  =  .6609  -  .4767  sin  (ws-1.047) 

and  w  is  the  hour  angle  and  ws  is  sunset  hour  angle.  In  equation  (A2) ,  R 
is  the  monthly  average  ratio  of  total  radiation  on  a  tilted  surface  to 
total  radiation  on  a  horizontal  surface;  therefore,  R  H  =  Hj,.  R  is  a 
function  of  both  the  beam  and  diffuse  radiation  incident  on  a  horizontal 
surface.  Since  these  measurements  are  generally  unavailable  they  are 
estimated  (as  discussed  in  Prins'  program)  at  hourly  intervals  from  H  by 
using  the  results  of  Liu  and  Jordan  (Ref  13).  They  found  that  the  ratio 
of  the  monthly  average  daily  total  radiation  (D/H)  was  related  to  K^,,  the 
monthly  average  long  term  clearness  index.  lCp  is  the  ratio  of  H  to  H  , 
is  the  extraterrestial  radiation  on  a  horizontal  surface.  To  find  the 
incident  beam  radiation  on  a  tilted  surface  when  the  incident  amount  is 
known  on  a  horizontal  surface,  the  ratio  R^  is  used.  R^  is  the  ratio  of 
beam  radiation  on  a  tilted  surface  to  that  on  a  horizontal  surface.  If  the 
tilted  surface  is  oriented  towards  the  equator,  then 


cos  (0-s)  cos^  cos  w  +  sin  (0-s)  sin*' 


cos  S  cos  w  +  sin 


where  w  is  the  hour  angle,  <f>  the  latitude  of  the  collector,  s  the  slope  of 
the  collector  from  the  horizontal,  and  S  the  declination  angle. 

The  utilizability,  IJT,  is  a  function  of  the  (to)  product,  the  average 
instantaneous  total  radiation,  I^t,  and  1J,  the  energy  loss  coefficient  of 
the  collector.  UT  is  the  fraction  of  incident  radiation  that  can  be  collected, 
or  utilized,  by  an  idealized  collector.  Utilizability  is  less  than  unity 
because  of  the  heat  loss  through  the  sides,  back  and  cover  plates  of  the  col- 
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lector  (Ref  16).  The  overall  energy  loss  coefficient,  U,  is  a  function  of 
the  collector  construction  and  its  operating  conditions. 

The  (7a)  product  is  a  function  of  the  dirt  factor  (DF) ,  the  shading 
factor  (SF) ,  and  the  angle  of  incidence  for  both  beam  and  diffuse  radiation. 

(tZ)  ■  .98  (1-DF) (1-SF) (to)  (A5) 

The  factor  of  .98  was  used  to  account  for  diffuse  radiation,  and  (xa)  is 
evaluated  for  beam  radiation. 

For  a  more  complete  and  detailed  discussion  of  energy  collection  one 
can  refer  directly  to  Captain  Prins  thesis  (Ref  16) . 
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Appendix  B 


Alternate  Methods  for  Determining  Heat  Transfer  Coefficients  and  Pressure  Prop 


Rock  Bed 

Lof  and  Hawley  (Ref  6)  also  investigated  packed  bed  energy  storage,  they 
arrived  at  the  following  expression  for  the  volumetic  heat  transfer  coeffi¬ 
cient  in  W/m^-°c,  hv: 


h 


v 


(Bl) 


2 

where  G  is  the  superficial  mass  velocity  in  K^/s-m  ,  and  D  is  the  equivalent 
spherical  diameter  of  the  particles  in  meters  given  by: 


net  volume  of  particles 
number  of  particles 


1/3 


(B2) 


As  before 


h  = 


h  AAx  h  D 

v  v  p 

A 


Another  empirical  correlation  for  packed  bed  heat  transfer  coefficients 


is  (Ref  4) : 

h  = 

and  jH  = 

where  the  Colburn  factor,  j^, 

E 

i 

i 

s 

and 


0.91Re'0,51; 

for 

Re  < 

50 

(B3) 

O.blRe"0,41; 

for 

Re  > 

50 

(B4) 

and  Reynolds  number,  Re,  are  defined  by: 

j„  =  CTJ  [Pr]2/*  (B5) 

P  ’ 
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Re  = 


lyV1^ 


(B6) 


where 

^8Ccn 

Pr  =  -  tt  *  =  Prandtl  number 

Kf 


G  =  pVQ  =  superficial  mass  velocity 

0  =  empirical  coefficient  which  depends  on  the  shape  of  the  particle 
(See  Table  Bl). 

A  6 

=  jp  (1-e  )  =  solid  particle  surface  area  per  unit  volume 
P  P  g 


The  subscript  of  f  denotes  properties  evaluated  at  the  film  temperatures ,  t£. 


Table  Bl:  Particle  Shape  Factors  for  Packed  Bed  Correlations  (Ref  4) 


Shape 

0 

Spheres 

1.00 

Cylinders 

0.91 

Flakes 

0.86 

Rasching  rings 

0.79 

Partition  rings 

0.67 

Berl  Saddles 

0.80 

Tube  Bundles 

A  tube  bundle  arrangement  may  also  be  used  for  thermal  storage.  The 
empirical  relationship  for  air  flowing  normal  to  a  bank  of  staggered  tubes 
is  (Ref  15) 

hi)  i  /, 

Nu  =  =  1.12  b7  (Rcj(Pr) 1  '  (B7) 

m  Rr  i 
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♦ 


where 


Re  =  Reynolds  number  =  pV  D  /ug 
1  max  o  - c 

V  =  velocity  based  on  minimum  flow  area 
max 

Pr  =  Prandtl  number  of  fluid 

h  =  mean  heat  transfer  coefficient 
m 

Dq  =  outside  tube  diameter 

Kj.  =  thermal  conductivity  of  air 


b^  and  n  are  the  Grimison  coefficients;  evaluated  from  tables  (Ref  15) 
and  dependent  on  tube  spacing. 

This  equation  normally  holds  true  for  tube  banks  with  ten  or  more  rows. 
However,  for  smaller  number  of  rows  a  correction  factor  exists;  these  factors 


Fig  Bl.  Staggered  Tube  Arrangement  (Top  View) 


Fig  B2.  Tube  Bundle  (Side  View) 
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IBS) 


V  =  Q/A  =  Q/D  L 
max  x  m  <  m 

where : 

A  =  D  L  =  minimum  flow  area 
m  m 

D  =  S™  -  D  =  minimum  clearance  distance 
m  T  o 

L  =  tube  length 

Q  =  volumetric  flow  rate 

To  determine  the  Grimison  coefficients,  b?  and  n,  two  more  parameters ,  Xj  and 

x™  must  be  defined. 

x,  =  Sj/I  )Q  (BO) 

and  x,j,  =  S.p/DQ  fB10) 


where : 


Sj  =  center  to  center  longitudinal  distance  in  direction  of  flow' 

S,f  =  center  to  center  transverse  distance  perpendicular  to  flow; 
coefficients  b^  and  n  arc  presented  in  tabulated  form  (Ref  15)  as  functions 
of  Xj  and  x,p. 

Equation  (B7)  is  valid  for  turbulent  flows  and  for  a  Reynolds  number 
range  of  2,000  to  40,000. 

A  pressure  drop  relationship  for  the  flow  of  air  over  a  bank  of  tubes 
also  exist  (Ref  8) : 


P  = 


f1G2max  N  „ 
j(2.09  x  Ifr ) 


(Bll) 


where 

G  =  V  =  mass  velocity  at  minimum  flow  area 
max  max 

n  =  density  evaluated  at  free  stream  conditions 
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N  =  number  of  transverse  rows 


v  =  wall-surface  temperature  dynamic  viscosity  of  air 
=  bulk-temperature  dynamic  viscosity  of  air 
and  the  empirical  friction  factor,  f\  for  staggered  tube  arrangements  is; 


f1  = 


0.25  + 


0.118 

Kst-d0)/i)0J^ 


-0.16 


Re 


max 


CR12) 


where 


pV  .  D 

Re  =  —  rnax- --  =  Refolds  number  based  on  minimum  flow  area, 
max  pg 
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Appendix  C 


Duct  Heat  Transfer  Relationships 

In  the  introduction  it  was  mentioned  that  the  inlet  air  temperature  for 
the  bed  during  its  charging  mode  was  approximately  the  same  as  the  temperature 
out  of  the  collector.  However,  it  is  only  when  duct  heat  losses  are  at  a 
minimum  can  this  assumption  be  made. 


Therefore,  to  validate  this  statement  a  prediction  of  heat  transfer  from 
circular  ducts  is  necessary.  The  following  diagram  is  used  for  this  purpose. 


where: 


Rj  =  forced  convective  resistance 

1*2  =  conductive  resistance  (insulation) 

=  free  convective  resistance 

For  the  fully  developed  turbulent  flow  region  the  forced  convective  heat 
transfer  within  the  circular  duct  can  be  represented  by  the  following  relation¬ 
ship  (Ref  10) . 


Nu  =  CRe 


.  08n  n 
Pr  = 


(Cl) 
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where  the  value  of  the  constant  is: 


C  =  .021  for  constant  axial  surface  temperature 
or  C  =  .022  for  constant  axial  heat  flux 
and 

h^  =  forced  convective  heat  transfer  coefficient 
=  thermal  conductivity  of  fluid  (air) 

Therefore, 

NUp  =  . 021Re‘ 8Prn  (C2) 

;  .4  <  n  <  .6 

Nil.  =  .  022Re‘  SPrn  (C.s) 


and  n  =  .6  is  commonly  used  for  circular  ducts  exposed  to  outside  ambient 
conditions.  Also,  these  two  algebraic  equations  are  valid  for  .5  <  Pr  *  1.0 
or  normally  for  gases,  such  as  air. 

The  first  two  resistances  arc: 


1 

R7 


R,= 


t. 

ins 

K. 

ins 


(C4) 


(C5) 


where 

t.  =  thickness  of  insulation 
ins 

K.  =  conduct ivi tv  of  insulation 
ins 

and  hj  is  determined  from  either  equation  ( (12 )  or  (C3).  Hie  recommended 
thickness  of  insulation  (Ref  1)  is  usually  one  inch  of  fiberglass  insulation 
for  duct  work  in  the  heated  area  and  about  1-2  inches  for  ducts  near  the 
collector  or  storage  unit. 
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To  complete  the  design  calculation  of  heat  losses  from  a  duct  to  an 
ambient  fluid  the  outside  resistance  due  to  free  convection  must  also  he 
taken  into  account.  For  free  convection  the  mean  Mussel t  number  can  be 
represented  by  the  following  equation  (Ref  8). 

Nun  «  C(GrDPrjm  (C6) 

where  the  constants  C  and  m  depend  on  geometry  and  whether  the  flow  is 
laminar  or  turbulent.  Also,  it  is  recommended  that  fluid  properties  should 
be  evaluated  at  a  mean- film  temperature,  t.. 

V  ■  *  tj/z 


where: 


t  =  wall  temperature 
t ^  =  ambient  (frec-stream)  temperature 
Also,  Gr^  =  Grashof  number  (based  on  diameter) 
where : 


gf 


(C7) 


and  >-■  =  coefficient  of  volumetric  expansion  =  r^-  ;  for  an  ideal  gas. 

1  op 


For  horizontal  ducts,  the  following  results  may  be  used: 


Gr„  Pr 


C 


m 


Approximate  Mean  Coefficients 


Laminar 


4  9 

10  -10 


.525 


1/4  ho  =  .27  $■ 


1/4 


CCS) 


Turbulent  10^-10^  .129  1/3 


ho  =  .18  (At) 


1/3 


(C9) 


where 


at  =  t2  -  t3 


D  -  “r  *  2  “ins5 

and  R_  =  —  ;  "ho"  can  be  found  from  either  equations  (CS)  and  ( C  9)  or 
3  Ko 

equation  (C6) .  The  total  heat  loss  relationship  can  now  be  expressed  as: 


QL  ~  UAs  (At)overall 

where:  U  =  overall  heat  transfer  coefficient  =  -  ^ 


(C10) 


or 


U 


*1 


+  Rt  +  R-j 


R 


eff 


and  A  =  total  surface  area  =  (Dl,  +  2  t.  3L...  Therefore,  the  heat  loss 
s  r  ins  D 

per  unit  length  can  be  written  as: 


QL  =  n(l>T  +  2  tins)  ^At)overal^ 

^  Reff 

where:  =  duct  length 

and  (At)  ^  =  tl  ’  t3  or  t^ie  overaH  temperature  difference. 

<J -  s  - * 


ceil) 


*id 


-fr 

t/qL 

Fig  C2.  Duct  Length  Section 
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cod 


can  also  he  expressed  in  the  following  form: 


ifi  c  (t.  ,-t  ,) 
a  pa1  id  od 


(Cl  2) 


where : 

=  air  temperature  into  duct 
tQ  j  =  air  temperature  out  of  duct 

Solving  equation  (C12)  for  "t  j"  we  get: 

t  i  =  t  ■  ,  -  Q.  /m  c 
od  id  'I-  a  pa 


(CIS) 


where  Qj  is  determined  from  equation  (Cl 1) . 

To  examine  the  changes  in  temperature  through  a  straight  duct,  several 
test  cases  with  variable  duct  lengths  were  carried  out.  In  all  cases  a 
negligible  temperature  difference,  tid'tod»  resulted.  Therefore,  the 
assumption  the  "t  j"  is  essentially  the  same  as  the  outlet  collector  tem¬ 
perature  proves  to  be  fairly  valid  for  our  particular  purposes. 
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Appendix  D 


Explicit  and  Modified- Implicit  Schemes 

The  governing  differential  air  equation  is  the  same  as  that  derived 
in  Chapter  II,  equation  (35). 


,T 
d  a 

dx 


(T  -T,  ) 
“acpa  a  b 


(M) 


Finite-differencing  equation  (Dl)  explicitly: 


T  ■  , -T  - 
a, l+l  a,i 

AX 


-KA 

ft  c 
a  pa 


(T  •  -T.  ) 
a,i  b' 


(h2) 


and  rearranging  gives 


T  h„AAxv. 

a,i+l  \  / 

’  \  a  pa/ 


a.i 


\,.\Ax 

ft  c 
a  pa 


C  H3) 


For  stability  purposes 


hvMx 

ft  c 
a  pa 


<  1 


or 


AX  * 


ft  c 


a  pa 

hvA 


p 


a^  ocpa 

'  hv 


(H4) 

(D5) 


Equation  (Dl)  in  a  modified- impl ic it  form  becomes 


Therefore 


T  •  ,-T  • 

a, l+l  a,i 

Ax 


-h„A 

ft  c 
a  pa 


-  T, 


(DO) 


a ,  i  + 1  - 


ivAAx  T 

ztnr  a, i 

a  paJ  ’ 


1  +  IkiAax 

Zfifc 
a  pa 
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WVax  ... 
ft  c  b 
+  JULa _ 

1  +  hyA^X 

2  ft  c 
a  pa 


(1)7) 


In  this  case,  stability  is  achieved  when 


IXjAax 

2m  c 
a  pa 


or 


Ax 


2m  c  2p  V  c 
4  3  Pa  _  a  o  pa 


hvA 


Here,  Ax  is  twice  that  of  the  explicit  result. 

For  these  two  schemes,  the  bed  equation  is  derived  in  the  same 
as  the  air  equation,  the  result  is: 


dT.  h  Aax 

“Ht  =  (1-Op  AAx~  (Ta'  V 

g  s 


(1-e  )p  c  A Ax 
g  s  ps 


In  explicit  form,  letting  B  =  rrr^O^c 

^  E  or  ^  ^ 


g  s  ps 


AAx 


equation  (DIO)  becomes 


T 


O  .  Tj 

b,i  lb,i 


nni,i-Tb,i>  - 


"L  At 


IT 


eS)pscPsMx 


or 


(1  -  B)  T 


BT- 


a,i 


0 


l  At 

(1-e  )p  C  .  AAx 
g  s  ps 


For  stability,  B  <  1  or 


At  <  n'e>scps 


Now  letting  C  = 


hyAAx 

me 
a  pa 


equation  (DA)  becomes 


(D8) 


(D9) 


manner 


(DIO) 


(Din 


(D12) 


(DIB) 
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gH9 

3‘t1 


Therefore,  equations  (D12)  and  (D14)  may  be  used  when  an  explicit  solution 
scheme  is  desired. 

The  bed  equation  in  a  modified- implicit  form  is 


t0+1 

b,i+l 


Tp 


b,i+i  - B  <V,i  -  rb,i+P 


'  At _ 

C1‘Eg)pscpsMx 


(D15) 


where 


*1 

a,i 


pj  +  t!  +  pj+l 

a,i  a,i+l  ^  a,i  a,i+l 
2  2 
2 


(D16) 


Substituting  for  T J .  from  equation  (D16)  into  equation  (D15)  and  rearranging 

a ,  i 

yields 


pj+l  _  pj 
b,i+l  o,i+l 


[1-B]  ♦ 


R/4  [T^  .  +  T-1  .  .  +  T^*  +  T^+I  ,1 
L  a,i  a, l+l  a,i  a,i+lJ 


(D17) 


Qi 


L  At _ 

(1-e  )p  C  Mx 
g  s  ps 


Since  C  =  ,  equation  (D7)  can  also  be  written  as: 


m  c 
a  pa 


Tj  -  U-c/2]  Tj  T  cl  i 

a,i+l  [i+t/Z]  a ,  i  |W2j  1b, 


i+1 


CD18) 


Therefore,  the  modified- implicit  solution  scheme  is  represented  by  equations 
(D17)  and  (D18). 
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Appendix  F 


Program  Instructions 

A  revision  of  Captain  Prins'  computer  program  for  solar  water  heaters 
was  necessary  to  deal  with  the  simulation  of  solar  air  systems.  The  modi¬ 
fied  program  can  handle  both  water  and  air  systems,  and  is  relatively  easy 
to  operate.  However,  when  an  air  system  simulation  is  desired  some  new 
data  cards  must  be  added;  consisting  of  rock  lied  and  pressure  drop  parameters. 

Whenever  the  air  system  is  in  effect  any  type  of  water  system  analysis 
is  overrided.  To  do  this,  A  PROG=4.0  data  card  is  all  that  is  needed;  the 
exact  location  of  this  card  will  be  further  explained  in  the  next  section. 

To  use  this  program,  one  must  also  input  data  about  the  weather,  posi¬ 
tion  and  orientation  of  the  collector  and  the  T.arth  as  explained  in  the  Prins' 
program  (Ref  16) . 

For  the  air  system,  rock  bed  parameters  are  used  to  determine  the  charac¬ 
teristic  heat  transfer  coefficients  between  the  air  and  the  rock;  the  pressure 
drop  parameters  are  used  to  determine  pressure  drop  and  blowing  costs  for  the 
bed,  collector  and  duct  work.  Next,  Prins'  part  of  the  program  is  used  to 
calculate  the  amount  of  useful  energy  collected  from  a  flat -plate  solar  collector. 
Once  this  is  done,  the  three  modal  operation  of  an  air  heating  system  can  lie 
analyzed. 

This  is  accomplished  by  calculating  a  space  heating  load  using  the  degree- 
day  model;  then  for  each  hour  of  the  day  the  useful  energy  collected  is  compared 
to  the  heating  load  requirement ,  and  based  on  this  criteria  a  particular  mode 
of  operation  is  chosen  in  order  to  satisfy  the  load. 

The  charging  and  discharging  modes,  modes  2  and  7>  respectively,  are  the 
most  important .  Here,  the  differential  equations  which  describe  the  storage 
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bed  are  needed  to  accurately  simulate  the  system.  As  explained  earlier  in 
Chapter  IV,  the  charging  mode  is  used  when  the  useful  energy  collected  is 
greater  than  the  space  load;  the  load  is  satisfied  first  and  anv  excess 
energy  is  stored  in  the  bed.  The  time  required  to  satisfy  the  load  and  the 
time  spent  charging  the  bed  arc  both  calculated.  Also,  if  the  discharging 
mode  was  previously  in  effect,  a  complete  bed  temperature  distribution 
reversal  is  performed  before  charging  the  unit. 

Mode  3  becomes  operational  whenever  solar  energy  is  no  longer  available 
for  collection,  in  this  case  energy  is  taken  from  storage  in  order  to  meet 
the  heating  load.  Again,  if  the  charging  mode  was  previously  in  effect  a 
bed  temperature  distribution  reversal  is  necessary  before  discharging.  As 
in  mode  2,  the  times  needed  to  satisfy  the  load  are  also  calculated. 

In  addition,  energy  losses  from  the  bed  are  calculated  for  each  period 
it  is  not  in  use;  updated  bed  temperatures  are  then  determined  at  each  nodal 
point  of  the  storage  unit. 

To  accurately  determine  fuel  savings,  a  yearly  space  heating  load  cal¬ 
culation  is  needed.  Also,  a  tax  credit  calculation  is  necessary  to  obtain 
the  net  investment  in  a  particular  type  of  solar  heating  system,  that  is 
total  investment  minus  the  tax  credit  break.  The  additional  parameters 
needed  for  this  cost  analysis  and  other  modifications  arc  fully  explained 
in  the  next  section. 

This  last  section  contains  the  extra  parameters  needed  to  operate  an 
air  system  simulation.  The  required  input  data  cards  arc  explained  along 
with  the  proper  nomenclature  and  input  format.  All  input  data  is  real 
except  where  noted.  Do  not  skip  any  data  cards.  Input  all  data  in  the 
following  order. 
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Card  # 

Symbol 

Input  Format 

Units 

Explanation 

1-15 

(First  15  data  cards  are 

the 

same  as  in 

Prins'  program,  Ref  16). 

16 

EFF 

x.xx 

none 

Efficiency  of  heat  exchanger. 
Air  systems  with  no  heat  ex¬ 
changer  enter  1. 

17 

PROG 

X. 

none 

If  PROG  =  4.0,  air  system  is 
in  effect.  For  water  system 
refer  to  Prins'  work  (Ref  16) 

18 

I  NT 

XX 

none 

Number  of  intervals  desired 
on  K~  curve  (INTEGER).  Ex¬ 
plained  in  Prins'  program 
(Ref  16) . 

19 

ND 

X 

none 

Number  of  consecutive  days 
the  simulation  is  to  run 
(INTEGER). 

20 

MON 

XX 

none 

Number  of  months  program  is 
run  (INTEGER). 

Cards 

21  through  46 

are  necessary 

for 

the  rock  bed  analysis. 

21 

TMIN 

XX. XXX 

minutes 

Time  increment  used  in  dis¬ 
charging  mode.  Used  in  a  Do- 
Loop  to  determine  when  load  i: 
satisfied  or  how  long  hot  air 
should  be  drawn  from  the  bed 
to  meet  the  load. 

22 

EPSIL 

XX.XXXXX 

hours 

Tolerance  use  to  leave  Do- 
Loop  after  1  hour  of  dis¬ 
charging  (1  hour  plus  or 
minus  tolerance) .  A  value 
of  .01  is  normally  used. 

23 

QMIN 

xxxx.xx 

B/hr 

Minimum  collected  energy  re¬ 
quired  before  fan  turned  on. 

24 

UA 

xxxx . xxx 

B/hr-F 

Overall  structure  conductance 

25 

DELIAS 

XX. XXX 

hours 

Desired  time  step. 

26 

Q 

XX. XXX 

cfs 

Volumetric  flowrate. 

27 

KB 

XX . XXX 

feet 

Bed  height. 

28 

WB 

XX. xxx 

feet 

Bed  width. 

29 

RHOA 

xxx. xxxx 

lbm/ft^ 

Air  density. 
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Card  # 

Symbol 

Input  Format 

Units 

Explanation 

30 

DP 

xx.xxxx 

lbm/ft3 

Rock  diameter. 

31 

PRN 

xx.xxxx 

none 

Prandtl  number  of  air. 

32 

CONDA 

xxx.xxxx 

B/hr-ft-F 

Conductivity  of  air. 

33 

NEUA 

xxx.xxxx 

ft2/hr 

Momentum  diffusivity  of  air. 

34 

CPA 

xx.xxxx 

B/lbm-F 

Specific  heat  of  air. 

35 

TOB 

XX.  XX 

Degrees -F 

Temperature  of  air  outside 
the  bed. 

36 

RHOS 

XXX. XXX 

lbm/ft3 

Density  of  rock. 

37 

CPS 

xxx.xxxx 

B/lbm-F 

Specific  heat  of  rock. 

38 

LB 

XX.  XX 

feet 

Bed  length. 

39 

KJ 

XXX 

none 

Number  of  nodes  in  bed 
(INTEGER). 

40 

EG 

XX. XXX 

none 

Void  fraction  of  bed. 

41 

KS 

XX. XXX 

B/hr-ft-F 

Thermal  conductivity  of  rock 

42 

PDIA 

XX. XXX 

inches 

Rock  diameter. 

43 

TBI 

XX.  XX 

Degrees -F 

Initial  bed  temperature. 

44 

KINS 

XXX. XXX 

B/hr-ft-F 

Thermal  conductivity  of 
insulation. 

45 

TINS 

XX. XXX 

feet 

Insulation  thickness. 

46 

AINS 

xxxx.xxx 

sq  ft 

Surface  area  of  bed  (top  and 
sides  only) . 

Cards  47  through  55 

are  needed  for  pressure  drop  calculations. 

47 

UG 

xx.xxxx 

lbm-hr 

Dynamic  viscosity  of  air. 

ftZ 

48 

LEFT 

x.xxxx 

feet 

Fan  efficiency. 

49 

ACC 

XX . xxxxx 

feet 

Collector  air  gap  spacing. 

50 

BC 

XX . XXX 

feet 

Width  of  collector. 

51 

LC 

XX.  XXX 

feet 

Length  of  collector. 

52 

LD 

XX . XXX 

feet 

Length  of  duct  work. 
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Card  # 

Symbol 

Input  format 

Units 

Explanation 

53 

DTD 

XX . XXX 

feet 

Duct  diameter. 

54 

TO 

XX.  XXX 

none 

Number  of  turns  or  elbows. 

55 

NH 

XX. XXX 

hours 

Hours  of  operation 

56 

DDAY 

XXX XX.  XX 

Degrees -F 

Number  of  Degree-Days,  (for 
the  average  number  of  Degree - 
Days  corresponding  to  each 
month  of  the  year  in  Vandalia, 
Ohio  refer  to  Table  £1). 

57-63 

(Same  as  cards 

28-35  in  Prins' 

program,  Ref  16) . 

Cards  64 

through  85  are 

for  the  cost  analysis  portion.  If  a  cost  analysis  is 

not  desired,  do  not  enter  any  more  data  cards  and 

automatically. 

the  program  will  stop 

64 

FTC1 

xx.xxxx 

none 

Primary  Federal  tax  credit 
percentage  (ex.,  40%  of  the 
first  2000  dollars  invested; 
per  cent  is  entered  in  decimal 
form,  40%  =  .40). 

65 

FfC2 

xx.xxxx 

none 

Secondary  Federal  tax  credit 
percentage  (per  cent  is  enterec 
in  decimal  form) . 

66 

OTC 

xx.xxxx 

none 

Ohio  tax  credit  percentage 
(per  cent  in  decimal  form) . 

67 

SICA 

xxxxx.xx 

dollars 

Initial  tax  credit  amount 
(i.e.,  first  2000  dollars). 

68 

LTCA 

xxxxx.xx 

dollars 

Secondary  tax  credi t  amount . 

69 

TMAYC 

xxxxx.xx 

dollars 

Maximum  return  from  tax  credit 
(i.e.,  amount  from  tax  credit 
can  not  exceed  4000  dollars). 

70-85  (Same  as  cards  36-51  in  Prins'  program.  Ref  16). 


Appendix  F 
Air  System  Flowchart 


Mode  Trigger 
Flow  2=0.0 
Flow  3=0,0 
JI=1 


PROG  FQ  4, 


READ  IN  BED 
PARAMETERS 


:ross  sec  are; 
Superficial  5 
Interstitial 
/elocities 


Mass  Flowrate 
Reynolds  H 
Heat  transfer 
coefficient 


Biot  # 

Effective  Coe: 
Constants  for  I  Cl  + 
Fin-Dif  Eqtnsf 


y _ 

Bed  Initial i 
for  each  pos 
tion  I)  time 


95 


UPDATE  OLD 
BED  TEM¬ 
PERATURES 


98 


Appendix  G 


Computer  Program  Listin' 
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